function getDifferentialOperators(self)
%Method to form differential & averaging operators
    
    n1 = self.ncx;
    n2 = self.ncy;
    n3 = self.ncz;
    d1 = ddx(self.hx);
    d2 = ddx(self.hy);
    d3 = ddx(self.hz);

    %Nodal derivatives
    D1  = kron(speye(n3+1),kron(speye(n2+1),d1));
    D2  = kron(speye(n3+1),kron(d2,speye(n1+1)));
    D3  = kron(d3,kron(speye(n2+1),speye(n1+1)));

    %Gradeint
    self.GRAD = [D1; D2; D3];
    
    %Cell center derivatives
    D1  = kron(speye(n3),kron(speye(n2),d1));
    D2  = kron(speye(n3),kron(d2,speye(n1)));
    D3  = kron(d3,kron(speye(n2),speye(n1)));

    % Divergence on faces
    self.DIV = [D1 D2 D3];
    
    %Curl on edges
    D32  = kron(d3,kron(speye(n2),speye(n1+1)));
    D23  = kron(speye(n3),kron(d2,speye(n1+1)));
    D31  = kron(d3,kron(speye(n2+1),speye(n1)));
    D13  = kron(speye(n3),kron(speye(n2+1),d1));
    D21  = kron(speye(n3+1),kron(d2,speye(n1)));
    D12  = kron(speye(n3+1),kron(speye(n2),d1));
    
    O1 = sparse(size(D32,1),size(D31,2));
    O2 = sparse(size(D31,1),size(D32,2));
    O3 = sparse(size(D21,1),size(D13,2));
    
    self.CURL = [O1     D32  -D23
                -D31   O2   D13
                D21   -D12   O3];

    %Averageing
    av = @(n) spdiags(0.5*ones(n+1,2),[-1,0],n+1,n);
    av1 = av(n1); 
    av1([1,end]) = 1;
    av2 = av(n2); 
    av2([1,end]) = 1;
    av3 = av(n3); 
    av3([1,end]) = 1;
    
    %Face avg operator
    self.Af = [kron(speye(n3),kron(speye(n2),av1)); ...
        kron(speye(n3),kron(av2,speye(n1))); ...
        kron(av3,kron(speye(n2),speye(n1)))];
    
    %Edge avg operator
    self.Ae = [kron(av3,kron(av2,speye(n1))); ...
        kron(av3,kron(speye(n2),av1)); ...
        kron(speye(n3),kron(av2,av1))];
    
    %Node avg operator
    self.An = kron(av3,kron(av2,av1));
end

function D = ddx(h)
    n = length(h);
    h = h(:);
  D = spdiags(1./[-h h],[0,1],n,n+1);
end